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Abstract 

Far-Held boundary conditions for external flow problems have been devel- 
oped based upon long- wave perturbations of linearized how equations about a 
steady state far held solution. The boundary improves convergence to steady 
state in single-grid temporal integration schemes using both regular-time-stepping 
and local-time-stepping. The far-held boundary may be near the trailing edge 
of the body which significantly reduces the number of grid points, and there- 
fore the computational time, in the numerical calculation. In addition the 
solution produced is smoother in the far-held than when using extrapolation 
conditions. The boundary condition maintains the convergence rate to steady 
state in schemes utilizing multigrid acceleration. 
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1 Introduction 

Numerical solution schemes for nonlinear flow equations usually require the introduc- 
tion of one or more artificial boundaries to be placed at some distance from a body 
around which (in external flows) the flow takes place. On these boundaries appropri- 
ate boundary conditions must be established that are not only physically correct, but 
that also comply with mathematical requirements as well. In, for example, subsonic 
flow regimes, it is well known that for hyperbolic flow equations (e.g. Euler Equa- 
tions), one characteristic points back into the flow domain at outflow, and therefore 
one condition must be prescribed; the rest, called numerical boundary conditions, 
must be in some way consistent with the partial differential equations. This paper 
deals with the derivation of outflow (far-held) boundary conditions. These boundary 
conditions will be both temporally and spatially local. 

The proper formulation of far-held conditions remains a vexing problem to date, 
and many authors have tried to tackle it in a variety of directions. It was once common 
in steady state calculations to set p = p^, (p being pressure) and to, say, extrapolate 
all other quantities. Certainly this procedure is inappropriate for boundaries close 
to the body, and so alternate approaches are needed to be developed. At present, 
the Navier-Stokes solver, [?], uses extrapolation on all physical quantities in subsonic 
regimes ( note that in supersonic regimes, extrapolation of all variables is correct from 
a mathematical standpoint ) with great success. In fact, its authors claim that they 
have been unable to replace far-held extrapolation with another far-held condition 
and to achieve satisfactory steady state convergence [?]. 

It has become popular to place “non-rehecting” boundary conditions at outhow 
in subsonic hows. While it is true that generally one cannot develop local boundary 
conditions that give no rehections, one can consider conditions which to some degree 
are better than others. The notion of better in this case implies that (see [?]) 

• the rehections are decreased and steady state convergence is accelerated. 

• the rehections will decrease as the position of the artihcial boundary tends to 
inhnity. 

• as the incident wave is more normal to the artihcial boundary, the rehections 
decrease. 

In 1988 Abarbanel, Bayliss and Lustman (ABL) [?] developed a non-rehecting 
outhow boundary condition for viscous, compressible hows over a hat plate. Their idea 
was to linearize the Navier-Stokes equations around a far-held steady state solution 
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Figure 1: External flow topology for compressible viscous flows around NAGA012 

with an assumed u-velocity profile based upon an asymptotic approximation to the 
ecpiations themselves. These linear perturbation ecpiations were then assumed to take 
on a dispersive wave solution and the resulting set of ecpiations (similar to the Orr- 
Sommerfeld system) described the asymptotic behaviour along the outflow boundary 
- the artificial boundary. Abarbanel, Bayliss and Lustman then assumed that the 
decay rate of the perturbations is mainly controlled by the iteration scheme’s inability 
to dissipate long wave disturbances. So they expanded the system, including the wave 
frecpiency, around small wave numbers i.e. long waves. The resulting zeroth and first 
order eigensystems are solved for their eigenvalues. The eigenvalue of the zeroth order 
eigensystem physically corresponds to the decay ratff: of the long wave disturbances, 
and the eigenvalue of the first order eigensystem physically corresponds to the group 
velocity of these disturbances. The slowest decay rate is found, and it defines unicpiely 
the group velocity. The details of the aforementioned strategy may be found in [?]. 
The outflow boundary condition derived from this takes on the form: 

ft + “(p - f-n = /k| cd 

where a > 0 is proportional to the decay rate of the disturbances, and (3 < 0 is their 
group velocity. This condition is very effective in the flat plate case. 

This paper uses the ABL approach in developing a far-field nonreflecting boundary 
condition for two-dimensional wake flows. We consider two cases. The first being 
compressible viscous laminar flows around NAGA0012, as illustrated in figure ??. In 
this case the compressible Navier- Stokes ecpiations are numerically integrated using 
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outflow 

boundary 


Figure 2: The flat plate topology for incompressible viscous boundary layer flows 

a time consistent scheme. Multigrid acceleration has been used in this case as well. 

The second case considered is incompressible viscous boundary layer flow over a 
finite flat plate as illustrated in figure •■'??. In this cas§ : the incompressible boundary 
layer equations are numerically integrated both time consistently and with local-tini© 
stepping. 

The nonreflecting boundary condition has been developed accordingly for time- 
consistent cases with and without multigrid acceleration and for local-time-stepping 
cases without multigrid acceleration. 

In §2 we outline., the long wave asymptotic expansion procedure in the general 
compressible case and develop the appropriate eigensystem for this case. In §3 w«. 
consider the incompressible boundary layer case and extend the boundary condition 
to local-time-stepping schemes. Then in §4 we derive the non-reflecting boundary 
condition both with and without multigrid acceleration. In §5 we present numerical 
results beginning with the numerical procedure used to obtain the eigenvalues needed 
in building the boundary condition, and ending off with a general discussion of the 
results of all the numerical tests performed. In §6 we comment on the results. 
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2 Long Wave Asymptotic Expansions 


The Navier-Stokes equations governing two-dimensional compressible flows may be 
written in the following form: 


continuity : 
x — momentum : 
y — momentum : 
energy : 


Pt + (pu) x + ( pv)y = 0 

p(,^t T au x T CUy^ T Px P-fl-xx T ^yy T x T ^y}x} 

p(v t + UV X + Wy) + Py = p{V X y + Vyy + | (u X + V y) y) 

P Cv A5t + P{ u x + v y ) = V ■ (kVT) + $ 


( 2 . 1 ) 


where p is the density, u and v are the x and y components of the velocity and p is 
the pressure. In addition, /c, the heat conductivity, c p and c V} the specific heats at 
constant pressure and volume, and their ratio, 7 = c v f c v are assumed to be constant. 
The viscous dissipation function, $ is defined as: 


$ — A (u x + v y ) 2 + 2 p{u 2 x + v 2 ) + p(u x + Vy ) 2 (2-2) 


and we shall use the Stokes relation, A + |/r = 0. 
gas, 


Cp c v p 


The equation of state for an ideal 


( 2 . 3 ) 


is used. We now define the Prandtl number, P r = — , which will be used in what 
follows. Also the Reynolds number with respect to the length of the airfoil is defined 
as 


Ri = 



( 2 . 4 ) 


where the subscript, 00 implies that the “free-stream” value of the quantity is taken. 
In addition, the viscosity, p is assumed to be constant. 


Assume that the flow is external, and the length of the airfoil is /. Defining the 
dimensionless variables: s = x/l and z = y //, an approximate wake profile in the zero 
lift case can be obtained [?]: 


u 


U n 


= 1 — m 


'/ri __ R ,£ 


( 2 . 5 ) 


Integrating (??) we obtain an expression for m in terms of the drag coefficient, 


CD- 


f +CO (i u w ° D 

m= (1 - 77 ) dz = W 


U a 


( 2 , 6 ) 
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We substitute the above into (??) and make the following parabolic approxima- 


tion: 


U 0 (z] S 0 ) def 
Loo = 


1 _ £D 

4V 7rs o 

1 


h _ 4:R,Z 2 1 -r 2 ^ 97TSn — 7 
/7TS0 97TSQ 4Rl P 


if 


( 2 , 7 ) 


where x 0 = s 0 l is far enough from the airfoil so that the value of to using (??) will 
not change as s > s 0 . Note that (??) satisfies (??), and defines an approximate wake 
thickness of: 

3 Firs 

6 = 2 \Ir,- 

Perturbing (??) in the region y 2 < S 2 around steady state conditions far down- 
stream gives: 


p't + u o(y)Px + PooK + PooV' = 0 


P<x>[ U 't + U o(y) U x + ( U 0 )y V> ] + p'x — p\ u 'xx + u 'yy] + f [ U 'x + V 'y ]a 
Poo[v' t + u 0 (y)v x \ + p' y = p[v' xx + v' yy \ + £[u' x + V' y ]y 


p't + 7Poo< + IPooVy + u 0 (y)p'x = 


IF 
Pr Pc< 


[{p'xx F p'y) — Fflip'xx + p'yy)\ + 2 (7 — l)(uo)yp[u' + V x ] 


( 2 . 8 ) 


where: v! = u — u 0 , p' = p — p 0} p' = p — p 0 and v' = u, and p 0 = Poo, Po = Poo, Ro = 0, 
u 0 (y) is dehned in (??). In addition, ^(-)o = §~ t (-)o = 0. 


Let 


K = 1 


Cdy/Ri 


(2.9) 


then 

uo(y) = U 00 [K + (l-K)P]. (2.10) 

Note that in the wake region, 0 < K < 1 and outside of the wake region, K = 1. 
Also, we make the following dispersive wave ansatz: 


IV 1 


FAv) 1 

u' 

e i(ipt+bx) 

Fi(y) 

v' 


Fslv) 

. p' . 


FAv) J 


Next, we non-dimensionalize in the following manner: 


F\ — PooG\ F 2 — LooGh 
f 4 = Poo uIg 4 F = u t L 

. V 


Fs — UooGo 
b = £ 


dy 6 dy 


(2.11) 


( 2 . 12 ) 
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where v = y/ Poo (kinematic viscosity). Substituting the ansatz and the non-dim- 
ensionalizations into (??) leads to the following set of dimensionless perturbations 
(similar to the Orr-Sommerfeld system) in the wake region ( — 1 < r] < 1). 


i[co + (5K( 1 + ( 1 r/ 2 )\Gi + iflG 2 + G' 3 — 0 

i\io + (3K(1 + ( 1 t/ 1 f] 2 )]G 2 + 2(1 — K)r]Gs + ifiG^ = s[G 2 + i 3 G 3 — pfd 2 G 2 ] 


i[co + (5K{ 1 + ? ? 2 )]G < 3 + — £ [|^3 + i 3 G' 2 — fd 2 G 2 ] 

i[to + /3/L(l + ( L ^ l V 2 )\G 4 + = 


(2.13) 


mG'i -$r)- 0*(G 4 - y%)] + 4(7 - 1)(1 - K)en[G> 2 + ifG z ) 


G 


where (•)' = 

We now perturb (??) around long wavelengths 
substitute: 

Gi = Gf ] + /3G\ 1] + !3 2 G (2) + • • • 

to = coq + /3wi + / 3 2 ^2 + ■ ■ ■ 


about small /3. To this end, we 


(* = 1,2, 3, 4) 


(2.14) 


into (??), obtaining the following zeroth and first order problems: 

The Zeroth Order Problem: 


' 9i- <p' = 0 

g 2 + 0g 2 — 2(1 — = 0 

<p" + !% - |V / = o 

g’l - 4 (-; - 1 )(1 - K)P,Ml<tA + W - tfUlW + “*] = 0 • 


( 2 . 15 ) 


where gi 



1,2, 3, 4 together with the change of variables: 


iu 0 = — eO </ ;•. = eQtp g\ = e 2 flt/j . 


( 2 . 16 ) 


The First Order Problem: 


Gi — — ^g 2 + ^[toi + K( 1 + lyy-?? 2 )]^ 

G" + ttG 2 - 20(1 - + A(1 + ^V)]</ 2 + *0e 2 [V> - §</] 

*" + - flT = if [a* + K( 1 + ^ 2 )]p - ^ 

G" - 4(7 - 1)(1 - K)P r MlgG' 2 + P r 0 $' - - \/ 2 0 II" + ^U] = 

-is 2 P r Ml 0[^ + K( 1 + 4^V)]<^ 

-iP r 92 + 4e 2 z(7 - 1)(1 - K)P r M^ o g0,<p 


(2.17) 


with the following change of variables: 


Gh = eGf } 


G* = eGP 


$ = ^3_ 

^ Q 


n = %- 


( 2 . 18 ) 
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In order to solve the zeroth and first order problems, appropriate boundary con- 
ditions need to be established. We shall derive these conditions by looking at the 
perturbations outside the wake region, K = 1, and then match at rj = ±1. One 
can easily check that placing K = 1 in (??) gives a system with constant coeffi- 
cients. We shall call the perturbations outside the wake region, Q, analogous to G 
in the wake region. Clearly this implies that Qi takes the form Qie rr] with Qi as 
constant. Now the exact boundary conditions for the zeroth and first order systems 
are Gf(±l) = QflH),i = 1,2, 3, 4. In order to get a relationship on the Gf’ s we 
use the fact that they satisfy ±1 - Since also Qi must satisfy 


[dv ]„=±i ~~ [ r( 5*]»=± i 

Gf(±l) are: 


i = 1,2, 3, 4, then clearly appropriate boundary conditions on 


dGi 

dg 


= \rGi 


77 = ±1 


\rj = ±l 


i = 1,2, 3, 4. 


(2.19) 


In order to obtain r in the previous expression, we do the following: 

Substituting K = 1 into (??) and using Qi = Qie rr] and hence Q\ = rQi 7 and 
Q'i = r2 Qi we obtain: 


+ (3)Q i + i(3Q2 + rQs — 0 
(i(u + ft) - er 2 )Q 2 - i/3^rQ 3 + if3Q 4 = 0 
(t(io + f3) - 43 er 2 )Q 3 - i/3^rQ 2 + rQ A = 0 
(i(lo + (3) - jtr£r 2 )Q 4 + p s M 2 r 2 Q 2 + ^rQ 2 + j^rQ?, = 0 . 

V i t j. y- iv j. cq iyj - oo 1J co 


We expand Qi and r around fd as in (??) and dehne 
zeroth order problem in the free stream region: 


qi and obtain the 


i(jj 0 q 1 + r 0 q 3 = 0 
(iujo ~ erl)q 2 = 0 
(uo 0 - ler 2 0 )q 3 + r 0 q 4 = 0 
. jrrlqi + r 0 q 3 + (uo 0 - rl^e)M^ > q 4 = 0 

where it is evident that either q 2 = 0 or ^ = — 0. We note that e rr] cannot 

decay as r] — > ±oo if = —0 hence q 2 = 0. For the rest of the system, the first, third 
and fourth equations in (??), we apply the existence criterion for nontrivial solutions 
of homogeneous systems, namely: 



/ 

ILOq 

O 

O 


det 


0 

i(jj 0 - I erl r 0 

= 0 


V 

e r 2 

Pr'O 

r 0 Mfliuo 0 - j^er 2 ) ] 



( 2 . 22 ) 
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We expand this around the first row retaining terms up to 0(e 2 Mfl). Assuming that 
I 7 O £ 2 Mfl <C 1 , then, we find that: 


r 


2 

0 


n 2 e 2 Ml 

m - 1 + (i - 


(2.23) 


For subsonic fiows with 7 = f.4 and P r = 0.72 and small e (physically reasonable) 
the latter result is negative. Hence the perturbations will not decay in 7 . Therefore, 
we choose: 

r 0 = TfLM^e (2.24) 

where the — sign is taken for positive 7 and the + sign is taken for negative 7 , so 
that the perturbations in free-stream will die out in 7 . 

The first order problem in the free stream region takes the form: 

iuoQ ^ + r 0 Q ^ = - i(ui + 1)71 - 7777 
(udo ~ erl)Q^ = *§*073 - *74 

< (zcu 0 - + roQi^ = (-*(^i + !) + |e*’o*’i)7 3 - ri7 4 (2.25) 

+ roQ^ + (iujo ~ j^erDM^Q^ = 

(-*(^1 + 1) + 7?erori)Af^74 - - *773 . 


The r 0 chosen above makes the coefficient determinant of the left-hand-side singular. 
Hence, solvability here requires that the same coefficient determinant with one of its 
columns replaced by the right-hand-side be singular as well. This criterion leads to: 


*7 — i*(^<7 + 1)4^00 + ()\ 2 M 2 !. 


(2.26) 


In order to fix the boundary conditions for (??) and (??), we recall the requirement 
that: 

= [rG t ] v=±1 * = 1 , 2 , 3, 4 

r; = ±l 

whose f3 expansion is 


dGi 

dr] 


d 

dr] 


[g«”> + edp + • • •] 

L J T] = ± 1 


7 + / Jr 1 + ...)(Gf> + / 3 Gf , + ---)' 

L J 7] = ±1 


gives 


/ Gf l '(±l) = r„Gf , 

t Gf»'(±l) = r 0 G,<"(±l) + r ,Gf> . (2 ' 27) 

This, together with 72 = 0 <-> 7 2 (il) = 0 provides us with appropriate boundary 
conditions for the zeroth and first order perturbation problems. The reader will note 
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that the systems (??) and (??) are of order six, so only six conditions will be 
prescribed for them. For physical reasons, we do not place conditions on the pressure 
perturbations, i.e on and on II. 

In summary, the boundary conditions on (??) are: 

g[(±l) = r 0 flfi(±l) 

<7 2 (±1) = 0 (2.28) 

<^'(±1) = r oF(±!) 

and on (??) are: 

Gy(±l) = r 0 Cri(±l) + £ri<7i(±l) 

G" 2 (±l) = r 0 G 2 (±l) (2.29) 

$'(±1) = r 0 $(±l) + eriLp(±.l) 

Since the perturbations decay like c ~ eQt , we solve the zeroth order problem for its 
minimal positive 0 and then use this 0 in the first order problem to find using a 
Fredholm like process. Since ^ |, 3 =o= <^i, then has the physical meaning of being 
the group velocity of the longest waves. 


3 The Incompressible Case 

In this section we shall find appropriate 0 and in the special case of incompressible 
viscous flows, i.e. M = 0. 


3.1 The Time Consistent Case 


The zeroth order problem in this case is: 

' 9i ~ ¥>' = 0 

g” + ttg 2 - 2(1 - K)riu<p = 

<p" + P<p - IV = 0 

g'{ + 0 P r Lp' = 0 

whose normalized solution, in the sense that || g 2 


n ^(±i) = o 
g 2 (± i) = o 
^'(±i) = o 

II r 2 =1 is as follows: 

11 (-M) 


gi = <p = i’ = o 


92(v) 


cos 



hence: 


0 < 0, 



(3.1) 


(3.2) 

(3.3) 
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The first order problem is interesting in that looking at (??) one would think that 
the appropriate boundary conditions on G 2 (g) are Gfdil) = 0 making the usage of 
the Fredholm alternative incorrect. It can be easily shown that in fact: 


Gi — — ^<72 

G" + n0 2 - 20(1 - K)r](f) 

$//+!„$_ * 11 ' = -^' 

, G'l + P r 0$' = -iP r g 2 


i[ u i + + L ^ L 'n 2 )\g 2 ^(±1) 

$'(±1) 


0 

0 

0 


(3.4) 


so that left-hand-null- vector of the first order problem is indeed cos(fy). Taking the 
scaler product of cos(fy) and the second equation in (??) gives 

(3.5) 

3 7w 

It should be made clear at this point that had the condition g 2 (i 1) = 0 been used 
instead of g 2 (i 1) = 0, in the zeroth order problem, we would have found that g 2 (g) = 
siring). In this case we would not have been able to solve the first order system for 
u>\. So it is very important to identify the proper boundary conditions in order to 
assure finding 0 and uq. 

The longwave expansion asymptotics done thus far are based on the assumption 
that the system of partial differential equations are solved in a time consistent manner 
using a global time marching scheme. In addition, had we started with the viscous 
incompressible boundary layer equations 


I u x + v y — 0 (3 6) 

\ U t + ( U 2 ) X + (uv)y = VUyy ^ ' 

and done the same type of longwave asymptotics as in §2 (note that the far-held 
profile, (??) is appropriate for this case as well) we would have obtained the same 
0 and loi as above. This means that essentially we have found the eigenvalues for 
the viscous incompressible boundary layer equations for time consistent numerical 
schemes. In the next subsection, we shall find these eigenvalues when local time 
stepping is to be used. 


3.2 The Local Time Stepping Case 

Many researchers use local time stepping in numerical marching schemes towards 
steady state. While local time stepping is not time consistent, it has been found to 
accelerate convergence to steady state. 

We return to the viscous incompressible hat-plate boundary layer equations, (??). 
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We are interested in integrating this system numerically. In order to ensure nu- 
merical stability, the Courant Fredrichs Levy (CFL) condition on the time step for 
each grid node requires (see [?]) that At be bounded by some function of A Xij and 
A yij. The function depends upon the numerical scheme used. In our numerical ex- 
periments Runge Kutta schemes were used. We were conservative in our approach 
and required that 

At < minA^-j (3-7) 

ij 


where 


A Oij < min 


o ( A S/«j) 2 

Uij + 1 ’ 2v 


( 3 . 8 ) 


where 0 < (3 < 1. The addition of the factor (3 helped stabilize the numerical iterative 
process. For the runs with global time stepping on a grid with constant Ax and A y, 
[3 = 1 was enough for stability, however for the runs on a grid with nonconstant A y, 
we chose (3 = | to achieve stability (see §5). 

Local time stepping means that at the node (x 8J , yij) a time step of A6ij is chosen. 
Using local time stepping one discretizes the temporal derivatives in the following 
manner: 



,n + 1 


~ U; 


AO; 


( 3 . 9 ) 


Suppose we discretize (??) on a non-uniform rectangular grid and temporally 
discretize as in (??). Provided that the spatial derivatives are discretized consistently, 
= 0, signifying that the grid is not changing in time. Therefore, 
these spatial discretizations are in fact valid approximations for the spatial derivatives. 
However, the temporal derivative is quite different. Since 


u l +1 ~ u [ 

A% 


' At ' 

AeTj, 


— vA- 


At 


we clearly see that what we are approximating is in fact 


lim 

At^O 


At 


du 

~di ' 


( 3 . 10 ) 


Suppose that the grid was constructed so that in the boundary layer, 

,(At/,) 2 


A % = (3 - 


2v 


(3.11) 


and that the grid is geometrically expanding in the y direction so that 


A y^ = k J Ayn] k > 1 
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then 

( A t \ = (A mV 

\A9ij J \Ayij J 


and therefore 


lim 


At 


At^o yA 0 {j 



= F(y,k). 


(3.12) 


Note that 0 < F(y, k) < 1 when k > 1 and that F(y, 1) = 1. 

Therefore, while we are solving (??) on this special grid using local time stepping, 
we are in actuality solving the modified system: 


U X + Vy = 0 

F(y; k)u t + ( u 2 ) x + (vu) y = vu yy 


(3.13) 


We shall now carry out a long-wave asymptotic analysis for the above system of 
equations. We note that the far-held wake profile, (??) is still appropriate in this 
case. Hence, using the parabolic far-held approximation for the wake prohle (??) we 
again dehne the primed variables 

u' = u — u 0 v' = v . 


Perturbing (??) around steady state variables retaining linear terms in primed vari- 
ables and substituting the dispersive wave ansatz: 


u' 

J(i pt+bx) 

' Fife) ' 

v' 

— c 

. Fl (s') . 


( 3 . 14 ) 


and nondimensionalizing gh , b } e, and y as in (??) we obtain the following system of 
ordinary differential equations in y: 


zGh(?7)/3 + G' 2 {y) = 0 

i[F{y k)u + /3K(1 + ^t? 2 )]^??) (3.15) 

+ 2G 2 (r,)(l - K)r, = eG’Kr,) 

where F\ = UooGi and F 2 = CooGh and 


F(ri;k) 



< 1 
> 1 


Perturbing Gi and u about f3: 

Gi = Gf ] + PGV + f3 2 G\ 2) + ■■■ (z = 1,2, 3, 4) 
to = COq T P^l T fFi 02 + • • • 


(3.16) 


(3.17) 
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and substituting these perturbations into (??) we get the following zeroth and first 
order problems in f3: 

The Zeroth Order Problem : 

v'iv) = ° 

g’l + 0 F(r\\ k)gi — 2(1 — K)rj£hp = 0 

where G\ 0 ^ = gi, g 2 = £^l<f and iaj 0 = — Oe. 

The First Order Problem : 

' ®' = ~hp 

< G’l + 0 F{rj- k)G 1 - 2fM(l - K)rj = 

m [F(r, k )^i + K{\ + (^) ri T )\ 

where G\ = eG ^ and T = 

Again, we need to match solutions at \g\ = 1 (at the boundary of the wake 
region) — when K = 1. Formally setting k = 1 because out of the wake region 
the grid stretching should have no effect either on 0 or on we obtain the linear 
perturbation system in the free stream region: 


(3.18) 


(3.19) 


f ! G 1 ( A3 + G)(9) = 0 
\ ![w + /J]G'i(!f) =eG’l( If) 

whose constant coefficients imply a solution of the form: 


G\ 

— ^r. 

' Gflr,) ' 

g 2 

— 0 

. g 2 (v) . 


Expanding as before, around f3 

Gi = gi + fdG] + ••• 
r = r 0 + f3r 1 + • • • 
lo = coo f3i0\ + • • • 


(3.20) 


(3.21) 


we obtain the following zeroth and first order problems in the free stream region: 


0 th order ((3°) : r 0 g 2 = 0 

(luo - erl)gi = 0 


(3.22) 
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from which immediately g 2 = 0. Since iio 0 = —eft then gi = 0 as well, otherwise the 
perturbations will not decay in g. 


1 st order (/3 1 ) : ig\ + r 0 G\ + rig 2 = 0 

G\(ito 0 - erl) = 0 

As in the zeroth order problem above, both G\ = 0 and G\ = 0. 


(3.23) 


This implies that the boundary conditions for the zeroth order problem are 


<7i(±l) = 0 

?(± 1 ) = 0 


(3.24) 


and that the boundary conditions for the first order problem are 

Gh(±l) = 0 

$(± 1 ) = 0 . 

Notice that the zeroth order problem reduces to: 

9" + fLF(g; k)g 1 = 0 
<7i(±l) = 0 


(3.25) 


(3.26) 


Since 0 < F(g;k) < 1 when k > 1 then the Sturm Comparison and Oscillation 
theorems (see [?]) indicate that the minimal positive eigenvalue ft in (??) will 
be greater than the minimal positive eigenvalue in the non-stretched incompressible 
case. This fact is crucial and to a certain extent indicates why local time stepping 
is so effective. Since it is believed that indeed the long wave perturbations are the 
slowest to converge to steady state, then we have just shown that the decay of these 
disturbances is faster in local time stepping regimes than in global time stepping 
regimes. 

In (§5) we present numerical evidence indicating that the nonreflecting boundary 
condition is quite effective in local time stepping regimes indicating that the absorbed 
long waves do in fact decay as we predict. This is an interesting illustration of one of 
the reasons why local time stepping is effective. 


4 Derivation of the Far-Field Non- Reflecting Bound- 
ary Condition 

Once the two eigenvalues, ft and are found, an appropriate far-held boundary 
condition can be developed to accommodate the physical situation at the far-held. 
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4.1 Derivation for Single-Grid Numerical Methods 

From the ansatz (??), we have: 

P = UlG, 

Substituting the dimensionless temporal and spatial variables: into the profile, we 
obtain (up to order fl)\ 

P = e^W+^+Wp^UlGtiri) + poo . 

Taking a derivative with respect to 6 and using the identity is 0 jto 0 = — yyO, we arrive 
at the far-held boundary condition: 


dp 4 , . dp 

_ = __ !)(p _ Pco)+ o 1 _. 

2 

In the incompressible case where 0 = ^-, the condition takes the form: 

dp 7 r dp 

^ = -g(?-Poo)+^. 


(4,1) 


(4,2) 


4.2 Derivation for Numerical Methods Utilizing Multigrid 
Acceleration 


The non-reflecting boundary condition just developed was based upon an analysis 
around long waves. Should multigrid acceleration be used in conjunction with the 
numerical scheme then short waves must be taken into account as well. Substituting 
P for p — poo, the boundary condition (??) may be rewritten as 


dP 

~dt 


4 

9tt 


91P + Ui 


dP 

dx 


(4.3) 


Taking the Fourier transform in x of the above equation we obtain 


Pt 


which has the general solution: 


P 


4 . 

— OP + iktoiP. 
9tt 


(4.4) 


For short waves, or high Fourier modes, k, the highly oscillatory behavour of the 
symbol P with t, will cause short wave disturbances in the numerical solution. This 
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is especially critical when using multigrid techniques whose goal on each grid is to 
damp short waves. Our original nonreflecting boundary condition, while effective in 
damping long wave disturbances caused by reflections, may create new short wave 
disturbances and destroy the effects of the multigrid acceleration. 


By replacing ik in (??) with 


ik 


1 + ik 

we obtain a new Fourier symbol free of oscillatory behaviour even as k gets large: 

P = . ( 4 . 5 ) 

The symbol (??) comes from the equation: 

P t = -—op + A Wl g . 

97 r 1 + ik 


Transforming back to the original variables we obtain the boundary condition to 
be used with multigrid acceleration: 


dp 40 dp d 

m = -^ {p - p ~ ) + u -''T x -d; 


% + h a{p ~ 


or 


d\ dp 40 

1+ aUS' = ~9T P ~ P “ ) + 




40 

9tt 


dp 

dx 


(4.6) 

(4.7) 


5 Numerical Results 

In this section we present various numerical results obtained using the far-held bound- 
ary condition just derived. We break this section into three subsections. The first 
describes the numerical technique for calculating 0 and in compressible cases. The 
second section deals with the numerical solution of the incompressible boundary layer 
equations both with global and local time stepping. The third presents the results 
obtained in compressible viscous flows around the NACA0012 airfoil using global time 
stepping regimes with and without multigrid acceleration. 


5.1 The numerical calculation of 0 and uj \ in the global time 
stepping case 

The discretized form of the zeroth order problem (??,??) is a generalized eigenvalue 
problem: 

Ax = h 2 klB(kl)x ( 5 - 1 ) 
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where h is the cell size, Ar/, and x is (g 1 , g 2 , ip, fl). The problem is non-linear since 
B = B( 0) is a function of 0, so the following iterative scheme was used: 

I 0 ° = eL 

{ Ax = h 2 tt n+1 B{tt n )x (5 ' 2) 

where at each iteration, the smallest positive 0 was determined. Generally conver- 
gence was obtained after three iterations and the main result is (Figure ??): 


For all physically relevant values of s, \l . and I\ ( see eq. ?? ); 

0 = f • 


The discrete first order problem (??, ??) can be written in the form: 

C(Q)X = (A- h 2 QB(Q))X = (5.3) 

where X = ( G\ , G-2 , $ , T ) and where A <->■ A and B B except that the boundary 
conditions on g- 2 are not parallel to those on &' 2 . The technique used to obtain u>\ 
was to find the left hand null vector of C Y (0), v, and then 


Ui = - 


( ^ ^i) 

(v, B 2 ) 


(5.4) 


Unlike 0, u>\ was not found to be constant for all physically relevant parameters. 
Figure ?? graphically shows the values of u>\ as a function of 1 — I\ for M iyo = 
0, 0.1, 0.2, 0.3 (for neatness \l . = 0.4, 0.5 were not shown). The main result is: 
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K-* 1 

uji — ->■ —1 for 0 < M < 0.5 and for all physical e. 

The limit \l . — > 0 is singular. 

In Table ?? values of 0 and u>\ for a variety of parameters, s, M iyo and I\ for 
time consistent numerical integration schemes are listed. These values were obtained 
through a program written in MATLAB. 


5.2 Incompressible Viscous Flow 

5.2.1 Case 1- Global Time Stepping Regimes 

The viscous incompressible boundary layer equations (for a flat plate) (??) have 
been solved in the geometry of Figure ??. The Reynolds number (with respect to 
the plate length), Ri, taken was 100000, and a rectangular grid (476x61) with cell 
aspect ratio, of 25 was used. The plate had 76 nodes in the x direction, and 
there were 400 x nodes in the wake region. At inflow, one quarter the way down the 
flat plate, a Blasius profile was placed. The outflow boundary was placed four plate- 
lengths from the trailing edge. The momentum equation was integrated in time using 
a Runga-Kutta scheme as suggested by [?], where the first order spatial derivatives 
were discretized using a compact fourth-order scheme developed by Abarbanel et. 
ah [?]. The second derivatives were discretized using standard second-order central 
differencing. 
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K 

e 

Moo 

0 


0.950000 

0.000100 

0.100000 

2.462296 

-0.979065 

0.950000 

0.001000 

0.100000 

2.462296 

-0.977492 

0.950000 

0.000100 

0.200000 

2.462190 

-0.979052 

0.950000 

0.001000 

0.200000 

2.462185 

-0.975585 

0.950000 

0.000100 

0.300000 

2.462014 

-0.978915 

0.950000 

0.001000 

0.300000 

2.461997 

-0.973317 

0.950000 

0.000100 

0.400000 

2.461773 

-0.978703 

0.950000 

0.001000 

0.400000 

2.461732 

-0.982999 

0.950000 

0.000100 

0.500000 

2.461471 

-0.977015 

0.950000 

0.001000 

0.500000 

2.461392 

-0.967744 

0.800000 

0.000100 

0.100000 

2.461776 

-0.911913 

0.800000 

0.001000 

0.100000 

2.461767 

-0.906617 

0.800000 

0.000100 

0.200000 

2.460254 

-0.904457 

0.800000 

0.001000 

0.200000 

2.460203 

-0.894001 

0.800000 

0.000100 

0.300000 

2.458007 

-0.897642 

0.800000 

0.001000 

0.300000 

2.457900 

-0.883191 

0.800000 

0.000100 

0.400000 

2.455157 

-0.892715 

0.800000 

0.001000 

0.400000 

2.455005 

-0.875650 

0.800000 

0.000100 

0.500000 

2.451718 

-0.892979 

0.800000 

0.001000 

0.500000 

2.451536 

-0.870731 

0.700000 

0.000100 

0.100000 

2.461120 

-0.862841 

0.700000 

0.001000 

0.100000 

2.461103 

-0.855154 

0.700000 

0.000100 

0.200000 

2.458011 

-0.847079 

0.700000 

0.001000 

0.200000 

2.457942 

-0.833530 

0.700000 

0.000100 

0.300000 

2.453518 

-0.840599 

0.700000 

0.001000 

0.300000 

2.453408 

-0.820463 

0.700000 

0.000100 

0.400000 

2.447673 

-0.835958 

0.700000 

0.001000 

0.400000 

2.447541 

-0.814361 

0.700000 

0.000100 

0.500000 

2.440383 

-0.834136 

0.700000 

0.001000 

0.500000 

2.440236 

-0.812478 


Table 1: 0 and u\ for various K, e, and M 
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Iterations 

Figure 5: Convergence behaviours for incompressible case. 


Two cases were run. The first with extrapolation of u at outflow, and the second 
using our far-held boundary condition (incompressible, (??)) at outflow with u>\ 
updated every iteration as in (??). In each case 5000 time steps were taken. Figure 
?? shows the convergence behaviour, ft is evident that the non-reflecting condition 
reduced markedly the residuals in the how held, implying that the the non-rehehecting 
boundary condition produced a more stable numerical solution than did extrapolation. 

In order to test the effect of the nonrehecting boundary condition as the outflow 
boundary is moved closer to the trailing edge, additional runs were performed with 
the outflow boundary at s = 2 and at s = 3. Figure ?? shows that as s decreases, 
the effects of the nonrehecting boundary condition are more dramatic. In addition, 
Call S2, S3, and S4 the solutions of the hatplate problem on grids G2 C G 3 C GA 
respectively. Use N as an extention for a case with the nonrehecting boundary condi- 
tions at outflow, and use E as an extention for a case with extrapolation conditions 
at outflow. We looked to see to what degree 

SN2 C SN3 C SNA, 


and 


SE2 C SE3 C SEA. 


It turns out that SN2 C SN3 C SNA up to four significant digits whereas SE2 C 
SE3 C SEA only up to two significant digits. Hence, we are certain that our bound- 
ary condition maintains accuracy better than the extrapolation condition does. All of 
this indicates that one could use a smaller computational domain while using the non- 
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centerline outflow 



centerline outflow 


Figure 6: Residuals after 5000 time steps - Extrapolation (top) Non-reflecting (bot- 
tom) 

reflecting boundary condition and maintain “large computational domain” accuracy, 
thereby dramatically reducing the amount of computational work needed. 

5.2.2 Case 2- Local Time Stepping Regimes 

We shall now describe the results obtained using the nonreflecting boundary condi- 
tion in conjunction with a local time stepping regime. Let us state the conditions of 
the case studied. 

(1) The Grid: As in subsubsection 5.2.1, the global time stepping case, the equa- 
tions are discretized on a 476 X 61 rectangular grid. The inflow boundary is located 
one quarter away down the flat plate and the outflow boundary is placed 4 chord 
lengths down stream from the trailing edge. Whiles' A.r is constant, Ay grows geomet- 
rically so that A yij = A;y 8j i, where A y hJ = y hJ — yij-i- The expansion factor, A’, 
was chosen so that in 60 intervals, y reaches one chord length. It should be noted that 
the aspect ratio of cells above the plate and the centerline are set identical to those 
in the global time step case so that in the local time stepping case At is determined 
by the parabolic part of the system in most of the wake region, recall (??). 

(2) Physical conditions: In this case, Ri is again taken to be 100000. Also, the 
boundary condition at inflow is the Blasius profile one-quarter way down the plate. 
The outflow boundary is located 4 chord-lengths away from the trailing edge. Equa- 
tions (??) and (??) with their boundary conditions, (??) and (??) are numerically 
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Ratio log(resid ext/resid nonreflect) after 5000 time steps 



t-,. « rrii , • 7 ( final residual with extrapolation bcA r , rnnn •. ,• 

rigure i: Lire ratio /ooi 0 E — i m — i — nn- 1 — i alter 5UUU iterations 

° J \ final residual with nrbc. ) 

as a function of the distance of the outflow boundary from the trailing edge of a flat 

plate. 


solved in this case for 0 and u>\ using a program written in Mathematica. The algo- 
rithm is similar to the the one presented in subsection 5.1. In this case, 0 is not a 
constant and depends upon the expansion factor, k. For these physical conditions we 
obtain 


0 = 5.95497 

2 

which is significantly larger than obtained in the global time stepping case. This 
implies that the long wave disturbances decay faster when local time steps are used 
than with global time steps. Regarding or, recall that in the global time stepping 
scheme case the relationship between between u>\ and I\ was linear. Therefore, the 
relationship between u>\ and q was linear as well. It has been numerically determined 
for the present local time stepping case, that the relationship between u>\ and q is 
linear as well, more specifically: 
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x 10 


-4 


(a) 


x 10 


(b) 


2 


4 



x 10 


-5 


(C) 


x 10 


(d) 


4 


4 



Figure 8: Residual plots after 5000 iterations (a) extrapolation, gts (b) non-reflecting 
condition, gts (c) extrapolation, Its (d) non-reflecting, Its 


or = 56.32 q - 2.01413 . (5.5) 

The non-reflecting boundary condition with these coefficients was used as an out- 
flow boundary condition, as was extrapolation. The results were compared and we 
state them now. 

In figure (??) the pointwise residuals after 5000 iterations for the four test 
cases (extrapolation outflow conditions-global time stepping, non-reflecting outflow 
conditions-global time stepping, extrapolation outflow conditions-local time step- 
ping, and non-reflecting outflow conditions-local time stepping), are plotted. In 
figure ?? the L 2 norm of the residuals (for each case) is plotted as a function of 
time steps. Figure ?? (a) graphs the logarithm of the L 2 norm of the residual as a 
function of time steps (iterations) for global time stepping runs. Figure ?? ( b ) graphs 
the logarithm of the L 2 norm of the residuals as a function of time steps for local time 
steping runs. The bottom dotted line is a special case that will be discussed later on 
in this section. It is clear from the results that 

• Local time stepping regimes converge faster than their global time stepping 
counterparts. 

• The case with non-reflecting outflow conditions with global time stepping reaches 
steady state faster than the case with extrapolation outflow conditions and local 
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(a) (b) 



iterations iterations 

Figure 9: Residuals (a) global time stepping (b) local time stepping 
time stepping. 

Note that adding the non-reflecting boundary condition to the local time stepping 
scheme speeds up the convergence. This implies that the non-reflecting boundary 
condition accelerates convergence to steady state above and beyond local time step- 
ping. 

For experimental purposes only we ran the flat plate code with local time stepping 
and the non-reflecting boundary condition using the appropriate 0 and u>\ for global 
time stepping schemes. The result is that the convergence behaviour in this case is 
unexpected. In fact, the convergence is better than with local time stepping with the 
non-reflecting condition along with it’s appropriate 0 and u>\. The bottom dotted line 
in (??) shows the convergence behaviour in this case, and another order of magnitude 
is obtained. Figure (??) shows the pointwise residuals after 5000 iterations in this 
additional case. Clearly th® best steady state is achieved in this way. 

This in no way disproves our claim. We obtain better convergence using our 
non-reflecting boundary condition than with extrapolation conditions both in global 
and local time stepping regimes. However, it appears that our method of hireling the 
coefficients is not optimal, so that perhaps factors other than long wave disturbances 
need to be taken into consideration. In any case, the important result that we have 
demonstrated here is as follows: 

Local time stepping schemes converge to steady state faster than global 
time stepping schemes do because they cause longwave disturbances from 
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Figure 10: Local time stepping run with 0 and lo\ taking on global time stepping 
values. Residuals 

steady state to decay faster than schemes utilizing global time stepping. 

5.3 Compressible Viscous Flow around airfoils 

5.3.1 Case 1- Global time stepping without multigrid acceleration 

The far-held non- reflecting boundary condition was implemented in the program [?], 
as a means to test its effectiveness in the compressible viscous subsonic case with 
Ri = 5000 and M iyo = 0.5. An appropriate 257 X 65 C-mesh was constructed using 
a hyperbolic grid generator. Often, hyperbolic generators create grids that tend to 
magnify the intrinsic singularity at the trailing edge. By adding additional control 
points we have overcome this phenomenon and generated a smooth grid reminiscent 
of an ellipticly generated one. The outflow boundary was placed 5 chord lengths 
away from the trailing edge, and was made ag perpendicular as possible to the airfoil 
chord. Two runs were made, one with extrapolation boundary conditions and the 
other with the nonreflecting boundary condition. In each run, 30000 time steps were 
taken. Neither multigrid nor other accelerators were used. Minimal artificial viscosity 
was used to assure convergence. The results that have been obtained are delineated 
below. 

The X 00 and L 2 norms of the p residual are shown in Figures ?? and ??. Clearly 
the nonreflecting condition accelerates convergence to steady state more than l lie 
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Residuals — Maximum Norm 



Figure 11: The X 00 norm of the p residual 

extrapolation condition does. The peaks in the residuals while using extrapolation 
have not yet been explained. Due to stability considerations (long time numerical 
integration), in both cases the residuals increase and level off towards the end of the 
runs, but in the nonreflecting case, the smallest residuals are smaller than in the 
extrapolation case, and the residuals grow sooner and faster with extrapolation than 
with the nonreflecting condition. Hence, the nonreflecting condition produces a more 
robust computational environment, and achieves a better steady state. 

In Figures ?? and ??, the pointwise residuals are shown for extrapolation and 
nonreflecting boundary conditions for 30000 time steps. While 30000 time steps is 
clearly beyond the optimal computation time, it is still evident that the nonreflecting 
boundary condition is more effective at reducing residuals than the extrapolation 
condition. In addition, in Figure ?? the values of p and u along the grid line: j = 5, 
i = 1 — 258 are shown. This grid line starts in the wake region (far held), then is in 
the boundary layer (around the airfoil), and returns to the far held. Both conditions 
give identical results around the airfoil - in the boundary layer. However, in the far 
held, the nonrehecting condition produces a smoother solution than extrapolating 
does. While we have shown this result for j = 5 only, it is true in the entire wake 
region. 


Far-Field Non-reflecting Boundary Condition 


August 2J h 1995 




Numerical Results 


27 


Average Residuals 



5.3.2 Case 2- Global time stepping withmultigrid acceleration 


In order to decrease the residual further and to stabilize the runs, we have investigated 
the nonreflecting boundary condition in regimes with global time steps, accelerated 
to steady state with multigrid techniques. 


As we previously mentioned, the nonreflecting boundary condition needs to be 
modified when using multigrid acceleration in the numerical scheme and it takes the 
form (??). 


In this case we compared two runs - one run using multigrid acceleration and ex- 
trapolation outflow conditions, and the other run using multigrid acceleration and the 
above nonreflecting boundary condition. In order to help eliminate short wave distur- 
bances in the y direction, we further introduced the Fourier smoothing as suggested 
by Ryaben’kii [?] 


P, 


hi 


:P\ 


* j+i 


H — P n ■ 

2 1 1 ’ 3 


+ i 


(5.6) 


locally near the lower and upper corners along the outflow boundary This smoothing 
transformation was applied each time our nonreflecting boundary condition was called 
and can be shown to identically kill off the shortest wave in the Fourier expansion. 


In the figures which follow (figures ??, ?? and ??), we show the convergence be- 
havior in both the L 00 and L 2 norms in using at outflow extrapolation and the nonre- 
flecting condition and show residual plots for p for both the extrapolation boundary 
condition and the nonreflecting boundary condition. 
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The same C-grid and other flow parameters were used in these calculations as 
in the previous subsection. Multigrid acceleration was applied in the following way: 
100 W cycles (3 grids) were run first starting on the first coarser grid in order to 
set up a reasonable initial condition. The results of these sweeps were used as the 
initial condition for a set of 2000 W cycles (4 grids) starting on the finest grid. Both 
second order and fourth order artificial viscosity (see [?] and [?]) were applied as 
well as implicit residual smoothing. Residual smoothing insures that the scheme will 
converge with a much larger GFL number hence accelerating convergence. 

The results are as follows: The L 2 residual norms obtained using both boundary 
conditions are nearly identical. However the residual plots indicate that the largest 
residuals at steady state with the both outflow conditions are at the corners of the 
outflow boundary (the corner problem) and at inflow, presumably propagated back- 
ward from the corners. The maximum norm, L 00 of the steady state residual while 
using the nonreflecting boundary condition at outflow, is slightly greater then when 
using the extrapolation conditions at outflow as is confirmed in figure ??. However, 
in the wake region, the residuals are smaller when using the nonreflecting boundary 
condition then when using extrapolation conditions. This is clearly seen in figures ?? 
and ??. 
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Figure 14: Residuals using Nonreflecting conditions for (a) p, (b) pu, (c) pv and (d) 

e 


6 Conclusion 

We have described the derivation of a nonreflecting far held boundary condition for 
steady viscous, compressible, two dimensional external flows. This condition pro- 
duces a smoother far held solution than the presently used extrapolation boundary 
conditions, while retaining the solution in the boundary layer. It also smooths out 
the trailing edge singularity by one to two orders of magnitude. The method can 
be applied to incompressible flows, where it is just as effective in driving residuals 
down to steady state faster than extrapolation conditions. In addition, nonreflecting 
boundary conditions of the type described in this work can be developed for schemes 
whose time stepping is local, as well as for schemes utilizing multigrid acceleration. 

In addition to the benefits accrued from using out boundary condition, it is easy 
to program, and does not require the use of much additional memory in an existing 
program. Also, its simplicity and local form require little cpu time. 

Perhaps the most important feature of our approach is that while in this paper we 
have chosen to deal with steady state problems, our methodology could be used to 
develop a condition for evolution equations that do not reach a steady state. In this 
case, we would perturb the linerized equations about the appropriate wave group that 
we wish to absorb (not necessarily about long waves) and proceed exactly as we did 
for steady state flows. Of course an appropriate far field profile needs to be formulated 
as well for each flow, which now may be time dependent. Hence our methodolgy is 
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Figure 15: The values of p and u along the grid line: j = 5, i = 1 — 258 


versatile and is shown to be effective in reducing residuals in external viscous flows. 
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log(average residual) 


log(max residual) 




Figure 16: The X 2 , (a), and L iyz \( b) norms of the p residual after 2000 multigrid W 
cycles. Legend: • • • - nonreflecting conditions, extrapolation conditions 
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Density Residuals- Muitigrid 



Figure 17: The p residuals after 2000 multigrid W cycles. Extrapolation outflow 
conditions. 


Density Residuals- Multigrid 



Figure 18: The p residuals after 2000 multigrid W cycles. Nonreflecting outflow 
conditions. 
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